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Non-axisymmetric oscillations of rapidly rotating relativistic stars are studied using the Cowling 
approximation. The oscillation spectra have been estimated by Fourier transforming the evolution 
equations describing the perturbations. This is the first study of its kind and provides information 
on the effect of fast rotation on the oscillation spectra while it offers the possibility in studying the 
complete problem by including spacetime perturbations. Our study includes both axisymmetric and 
non-axisymmetric perturbations and provides limits for the onset of the secular bar mode rotational 
instability. We also present approximate formulae for the dependence of the oscillation spectrum 
from rotation. The results suggest that it is possible to extract the relativistic star's parameters 
from the observed gravitational wave spectrum. 

^ PACS numbers: 04.30.Db, 04.40.Dg, 95.30.Sf, 97.10.Sj 
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I. INTRODUCTION 



Relativistic stars during their evolution may undergo oscillations which can become unstable under certain condi- 
tions. Newly born neutron stars are expected to oscillate wildly during their creation shortly after the supernovae 
collapse [T] . They are expected also to oscillate if they are members of binary systems and there is tidal interaction 
(3JT) [5J [3] or mass and angular momentum transfer from a companion star and also when they undergo phase transitions 
[HE] which might be responsible for the observed glitches in pulsars. Rotation strongly affects these oscillations and 
perturbed stars can become unstable if they rotate faster than some critical velocities. During these oscillatory phases 
of their lives compact stars emit copious amounts of gravitational waves which together with viscosity tend to suppress 
the oscillations. The oscillations are divided into different families according to the restoring force [BJ|7]. If pressure 
is the main restoring force then these modes are called (pressure) p-modes, buoyancy is the restoring force of another 
^sO class of oscillation mode, the g-modes while Coriolis force is the restoring force for the inertial modes. Spacetime 
induces another family of oscillations which couple only weakly to the fluid, these are the so-called w-modes |8J. There 
0^ are more families of modes if one assumes the presence of crust [BJ |SJ EH HI] or magnetic fields [T5] . For a complete 
description of the relativistic star perturbation theory one may refer to [7] [13j [14] . 

The study of stellar perturbations in the framework of general relativity (GR) dates back to mid 1960s with the 
seminal works by Thorne and his collaborators (HO EE] EO EE] . Since then the study of stellar oscillations and possible 
;> instabilities was a field of intensive work in relativistic astrophysics [TO] US]- Moreover, during the last two decades, 
these studies become even more important due to the relations of the oscillations and instabilities to the emission of 
gravitational waves and the possibility of getting information about the stellar parameters (mass, radius, equation 
of state) by the proper analysis of the oscillation spectrum [SU] [5TJ [55] [531 US (H] [5B]. Still, all these studies were 
mainly dealing with non-rotating stars, because the combination of rotation and general relativity made both the 
analytic and numerical studies extremely involved. This led to certain approximations in studying rotating stars in 
GR. The most obvious of them include the so called "slow-rotation" and the "Cowling" approximation. Actually, 
both approximation were known and have been used extensively in Newtonian theory of stellar oscillations |27j . In the 
slow rotation approximation one expands the perturbation equations in terms of a small parameter e = £1 /£Ik , where 
fl is the angular velocity of the star and £Ik stands for the "Kepler angular velocity" which is the maximum velocity 
that can be attained before the star splits apart due to rotation. In the Cowling approximation one typically neglects 
the perturbations of the Newtonian potential or the spacetime in the case of GR. This is a quite good approximation, 
both qualitatively and quantitatively, for the higher order p-modes, for the g-modes and the inertial modes while it is 
only qualitatively good for the f -modes, see for example |28j . Although, most of our understanding on the oscillations 
of relativistic stars is due to perturbative studies, recently, it has become possible to study stellar oscillations using 
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evolutions of the non- linear equations of motion for the fluid [551 USB [3D E21 1331 13H EUD [5SJ [37] • Finally, differential 
rotation is another key issue that is believed to play an important role in the dynamics of nascent neutron stars. 
Actually, it is associated with dynamical instabilities both for fast and slowly rotating neutron stars [3H1 [3H] and 
affects the onset of secular instabilities; still it has not yet been studied extensively [101 SU H2] and remains an open 
issue. 

Since rotational instabilities are typically connected with fast rotating stars, it is of the great importance the study 
of oscillations of this type of objects. As a first step one can drop the "slow rotation" approximation but still maintain 
the Cowling approximation i.e. to freeze the spacetime perturbations or in the best case to freeze the radiative part of 
these perturbations. This was the approach used up to now for most of the studies of the oscillations of fast rotating 
relativistic stars either in perturbative approaches [43j[44] or in non-linear (but axisymmetric) approaches [35l 145], 

In this article, we present the first results of a new approach based on 2D time evolution of the perturbation 
equations which seems to be promising for the study of the oscillations and instabilities of fast rotating neutron stars. 
This the first study of its kind, since earlier 2D perturbative studies have been done either in Newtonian theory 
[46] (see also [47] for a very recent work) or in reduction to eigenvalue problem [43l HH [48]. The advantage of this 
method is that it can be easily extended to include differential rotation or the perturbations of the spacetime. On 
the other hand it provides a robust tool in studying the onset of rotational instabilities in fast rotating neutron 
stars, while as it has been demonstrated here, one can get easily results for realist equations of state which is vital 
for developing gravitational wave asteroseismology. Finally, via this approach one may answer the question of the 
existence or absence of a continuous spectrum for the inertial modes as it has been suggested by the ID studies in 
the slow rotation approximation SHI (SOI [511 • 

In the next section we describe in detail the derivation of the perturbation equations and the conventions that we 
have adopted. In section 3, we present the numerical tools that have been developed in order to study the problem 
while section 4 describes the results for axisymmetric and non-axisymmetric perturbations. In the last section we 
discuss the results, their application to astrophysics and the possible extensions of this work. 



II. THE PERTURBATION EQUATIONS 

The study of oscillations of rotating neutron stars involves the solution of the full nonlinear set of Einstein's 
equations of General Relativity together with the equations of motion for the fluid (we set G = c = 1 here) 

G^v — (1) 



V^T"" = 0, (2) 

where G M „ is the Einstein-tensor, describing the geometry of space-time, and is the energy-momentum-tensor 
that defines the functional form of energy momentum and stress of the fluid. Since it is very complicated, but not 
impossible [35 , 45j, to solve this system as such, we will introduce some approximations in studying the problem. 

First of all, we linearize equations ([lj and which means that we constrain our study to small perturbations 
around the equilibrium. Secondly, we will work in the so-called Cowling-approximation, which means that we will 
neglect all metric perturbations. This simplifies significantly the equations one has to solve since the space-time is 
considered as "frozen"; and we only have to solve the linearized version of pi). Under this assumption equation ^ 
will be written as 

V^TH = g^{d^6T KV - T x K ^5T Xv - T x ^5T kX ) = (3) 

where g^ v is the metric describing neutron star's spacetime, T x KfJ _ are the Christoffel symbols and ST^ is the Eulerian 
perturbation of the energy- momentum-tensor. We assume that the matter has no viscosity or shear stresses, i.e. that 
it can be described by a perfect fluid. Thus ST^ has the form 

ST^ = (e +p)(u ll 5u u + u v 8u^) + (Sp + 5e)u^ + bpg^ . (4) 

where e is the energy-density, p is the pressure, is the 4- velocity and 8u^ are its perturbation. Energy density and 
pressure are not independent quantities but are connected via an equation of state (EOS) which we assume to be 
polytropic, i.e. 

p = Kp 1+1/N where e = p + Np. (5) 

Here p is the rest-mass density, K the polytropic constant, N the polytropic exponent and T = l + l/N the polytropic 
index. For barotropic oscillations, both the unperturbed background and the perturbations are described by the same 
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equation of state. In this case the pressure variation can be replaced by the corresponding energy-density variation 
via Sp = c 2 Se, where c s is the speed of sound 



C S =g~ e . (6) 



For polytropic EOS of the form ([5| it is given by 



(7) 



P 



Our background model is a compact relativistic star that rotates uniformly up to its Kepler-limit, i.e. the point were 
it is torn apart by centrifugal forces. In this work we will adopt the metric in a comoving frame of reference as it is 
described in |53j . In Lewis-Papapetrou coordinates (p, <p,t) the metric reads 

ds 2 = e- 2U [e 2k (dp 2 + d( 2 ) + W 2 dp 2 ] - e 2U (dt + adp) 2 , (8) 

where all functions depend on p and £. In general, the (ip, i)-metric component is proportional to the function a and 
vanishes in the absence of rotation. Other properties of the metric functions which will become important later are 

lim |ap~ 2 | < oo and lim | Wp~ 1 1 < oo . (9) 

The utilization of a comoving frame of reference as implemented in |53j has several important advantages. First of 
all is the use of surface-fitted coordinates. Especially when the star is rapidly rotating, deviations from spherical 
symmetry become more and more evident. If one uses a fixed spherical symmetric grid to describe the neutron star, 
even if the boundary of the nonrotating configuration can be aligned to a single grid line, this will fail once the 
star rotates and gets deformed, so that the surface will lie somewhere between different grid lines. This may cause 
problems when implementing boundary conditions at the surface. With surface-fitting coordinates, even the surface 
of the most rapidly rotating configuration can be described by a single parameter. Secondly, comoving coordinates 
reduce the complexity and length of the equations to solve considerably. Let us briefly review how one would proceed 
in a stationary frame: The first observation is, that there are still dependent quantities in the equations ^ and Q. 
One would write down the definition of the constant angular velocity 

u v 

— = and use the relationship Qu^u^u" = — 1 
u l 

to get u* = (— gu — 2f2<7 v t — f2 2 g w ) -1 / 2 . This is how the angular velocity enters the equations explicitly. Furthermore 
since also 



gtiuiu 1 * + Su^u" + Su") = -1, we have Su t = -Q5 



In this way we can write both it* and <5u* as functions of the angular velocity. On the other hand, in a comoving frame 
u v = while the equation for 5ut is trivial; i.e. Su t — 0. This reduces the equations governing linear perturbations 
considerably as we will see now. 

For non-rotating or slowly rotating stars where the background configuration was considered as spherical the 
perturbation equations were decomposed into spherical harmonics and the problem was typically reduced in solving 
the equations describing only the radial components of the perturbations. But here we deal with fast rotating neutron 
stars which are deformed due to rotation, this means that it is no longer possible to decompose the angular part of 
our perturbation quantitites in spherical harmonics as it was usually done for the non-rotating case (see [15 ) or in 
the slow-rotation-approximation (for example [54 and [55 ). Instead we can only separate the azimuthal dependence 
and the perturbation functions that we use will be written as 

(e + p)We u Su p = h(p,<:,t)e lmv 

(e + p)We u Su c - f 2 ( P X,ty m ^ (10) 
(e + p)Su v = MpX,t)e lm ^ 
c 2 e u Se = H(pX,ty m *. 

The functions fi ,i — 1 ... 3 are describing the time-evolution of the perturbed velocity components and H describes 
the corresponding change in energy-density. All these functions are in general complex-valued (of course not for 
the initial starting condition at t = to) and multiplied by a complex phase that prescribes the dependence on the 
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azimuthal angle ip. This means, that we will get a set of complex-valued PDEs and since this is a linear system, the 
final solution is obtained by taking the real part of these quantitities. 

The substitution of the perturbation functions (10 1 into equations ([3]), leads to the following system of evolution 
equations: 



where 
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(ii) 



F 



a 2 c 2 e 4U 



W 



(12) 



As we discussed earlier, there is no explicit dependence on the angular velocity £1 in this system of equations. 



The perturbation equations (11) are complemented by boundary conditions which describe the behaviour of the 
perturbations on the boundaries of the numerical domain, which are the rotation axis and the surface of the star. 
One also has to discriminate between perturbation variables with scalar and vectorial character to find the correct 
conditions along the rotation axis. In our case, the energy density perturbation which is described by H and the 
^-component of the perturbed 4- velocity Su^ 7 described by h are clearly the scalar perturbations while j\ and h 
(describing velocity perturbations in a ( = consi.-plane) are the vectorial perturbations. Let us first look at the 
boundary conditions at the surface, since this is very easy in our formulation. Since the perturbation functions on 
the lcfthand side of (10 1 drop to zero there, all our perturbation quantities vanish. So on the surface of the neutron 
star we have 



/ 



1 1 surface 



surface 



surface 



= H 



surface 



= o. 



(13) 



For the rotation axis we have to consider three different cases, depending on the value of m and the nature of the 
perturbation variable. Scalar perturbations have to be unique along the axis for all values of m when varying the 
azimuthal angle ip. Vectorial perturbations have to change like cos(tp) or sin(y>) when varying ip. This means, that 
only for m = ±1 they are allowed to have nonzero values along the axis. If we then take into account again the 
equations ^ and (10), we end up with the boundary conditions for the rotational axis depicted in table [I] 



m-value 


fl axis 


fl | axis 


^3 | axis 


H | axis 














finite & continuous 
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finite & continuous 





else 















TABLE I: Boundary condition for the perturbation variables along the rotational axis 



III. NUMERICAL METHOD 



As described already in the previous section, in this study we adopted a comoving frame of reference, in which the 
metric takes the form shown in equation ^ . The numerical method used to solve the equations governing the stellar 
background for this special metric is described in detail in |53j ; here we will briefly summarize the parts which are 
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crucial for our work. Since the stationary background model possesses rotational symmetry as well as symmetry with 
respect to the equatorial plane, for the computation it is sufficient to consider the physical domain 



z>+ = [(p,C),p>o,c>o] 



By means of a coordinate transformation 



[(a, r) G [0, 1] x [0, 1] , (p = p(a, r), ( = ((a, r)) G V+) 



(14) 



(15) 



the neutron star interior is mapped onto the unit square and then, the equations are solved in this new coordinate 
system with a spectral-methods-code. 

For example, the distribution of the grid points is shown in Figure [l] for a resolution of 18 x 18. Note, that the 
special choice of collocation points guarantees that there will never be any grid points directly on the boundaries. 
One can also observe the typical clustering of the grid points near the boundary which is characteristic for spectral 
methods. Also in this figure, some properties of the coordinate transformation T + are labelled. One can notice, that 
the surface of the neutron star gets mapped onto a — 1; this is independent of the rotation rate, even for rapidly 
rotating neutron stars this coordinate line always corresponds to the stellar surface. Additionally, r = corresponds 
to the rotation axis above the equatorial plane and t — 1 to the equatorial plane itself. As mentioned earlier, since the 
background star is axisymmetric as well as symmetric with respect to the equatorial plane, this computational domain 
suffices in order to construct the background model. Note also, since we have four boundaries in our computational 
domain, but only three "physical" boundaries (i.e. rotation axis, equatorial plane and stellar surface), one single point 
gets smeared into a coordinate line; in our case, a — corresponds to the center of the star (i.e. p — C, = 0). The 
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FIG. 1: Layout of the numerical grid used in constructing the stationary neutron star background model 



coordinates (a, t) are similar to their spherical coordinates counterparts (r, 9) in the sense that moving from r = to 
t = 1 on an arbitrary a = const-line means to start from a point at the rotation axis and move somewhat "parallel" 
to the surface to a point at the equatorial plane. Vice versa if one moves along an arbitrary r = const. -line then one 
starts from the center of the star to the surface. Figure [2] illustrates this for a rapidly rotating stellar model with a 
ratio of polar coordinate radius to equatorial coordinate radius of r p /r e = 0.6. This compact star is rotating near 
its Kepler-limit and obviously is non-spherical. Nevertheless, as discussed above, the surface of the stellar model is 
described by S = [(p(l, r), ((1, r)) , r G [0,1]]. 

For the linear perturbations under study these coordinates suffice and we used it with minor modifications. Because 



of the equatorial symmetry of the background it is sufficient to use T> + (see ( 14 1 ) as the domain of computation, but 
this is no longer true for general axisymmetric perturbations. There are two ways to circumvent this. The first is 
based on the fact, that every perturbation can be decomposed into a symmetric and into an anti-symmetric part. 
To illustrate this, consider A.+ — [x , x > 0] as a one-dimensional computational domain and f(x, t) as an arbitrary 
perturbation with values in A.— — [x , x < 0] too. If x = corresponds to the axis of symmetry we have 



with 



f(x,t) = f s (x,t) + f a (x,t) 



f.(x,t)=H X > t) + H- X > t) and /. M = M^ 
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FIG. 2: Some r = const, (solid) and a = const, (dotted) coordinate lines in the {p, £) coordinate system 

The procedure for a time-evolution of an arbitrary initial perturbation f(x,to) on A := A.+ + A- would be to 
decompose any initial perturbation into its symmetric and anti-symmetric parts, perform an evolution of these parts, 
complement the solutions according to their symmetry behaviour for x < and add the solutions up. Note that in 
this case there are two different boundary conditions to impose on the unknown function at x = 0. For the symmetric 
part we have d x f s (x, t)\ x =o = 0, for the antisymmetric part we have f s {x,t)\ x= o — 0. Note also, that in order to 
study eigenmodes and compute their eigenfrequencies it is sufficient to restrict these studies on either symmetric or 
anti-symmetric perturbations; for this purpose it is not neccessary to consider arbitrary initial data. 

Since the goal was to create a code that can handle arbitrary perturbations, we chose the second option in handling 
this problem, that is we extend our physical domain to include also the region 

2?- = [(p,0 > /»>0,C<Q]- (16) 

In practice we "glue" two copies of T> + together along the equatorial plane. Correspondingly the computational 
domain extents in the r-dimension beyond r = 1, i.e. the ^-direction in analogy to spherical coordinates we discussed 
earlier. Hence the extended domain is given by 

T = [(a, t) e [0, 1] x [0, 2] , (p = p(a, r), C = C(*> r)) £»:= (V+ + 2?_)] (17) 

With this choice, there are no other boundary conditions to impose than those described in section [TTJ the equatorial 
plane now lies in the interior of the computational domain and its boundaries are given by the center (<r = 0), the 
surface (a = 1), the part of the rotation axis above the equatorial plane (r = 0) and the corresponding part underneath 
the equatorial plane (r = 2). Of course one pays a price for studying arbitrary perturbations: This grid is now twice 
as large as before; this means for a 2D-code an increase in computation time by a factor of 4. 

Special care has to be taken for the correct transformation behaviour of our various background quantitites which, 
up to now, are only known in T> + . In addition to simply "mirror" all neccessary quantities along the equatorial plane 
(i.e. along r = 1), some derivatives need to be multiplied by a factor of —1. Derivatives with respect to p pose no 
problem; if (a, t) denotes a grid point in T> + (i.e. r < 1) and g is an arbitrary background quantity, we have 



9g, _ dg 

dp^ T) ~ dp 



(<r.T) - -^r\(a,2-r) ■ (18) 



This is no longer true when one uses derivatives with respect to £; in this case it is 



Qj\{a,r) - —Qf\(v,2-T) ■ (19) 



Yet another transformation behaviour has to do with the fact, that the equations we study (i.e. system (111) are 
written in (p, (^-coordinates but we use as numerical domain our extended (<r, r)-system. For our perturbation 
variables we have to switch between these systems. So in addition to all background quantities and their derivatives 
with respect to p and £, the transformations coefficients da /dp, da/dQ, dr/dp and dr/dC, are also available in T> + . 
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With these coefficients it is possible to compute for any perturbation quantitiy / given on the numerical domain T 
the values of 

(20) 



df = 


df da 


df dr 


dp 


da dp 


r — 

dr dp 


df 


df da 


dfdr 


d( 


dadC ' 


~ &rd( 



we need to know for the righthand-sides of the evolution equations (11). Equations ( 20 I are obviously valid for the 



background quantities as well and this helps in finding the correct form of the transformation coefficients when going 



from T> + to X>_. We know how the p- and ^-derivatives (i.e. the lefthand-sides of (20 1 ) transform in X>_ as well what 
happens to their a- and r-derivatives there. This leads to 

da da da da 

dp' (<7,T) = ~dp~* {a ' 2 ~ T) g^k'r.T) = --^:|(<T,2-T) (21) 

and 

dr. dr dr dr 

dp~k<r,T) = -q-\(<t,2-t) and q£\(v,t) = Q£k<r,2-r) ( 22 ) 



We then use standard finite-differencing schemes and time integrators to solve system (111. However, the numerical 
evolution of the equations was unstable i.e. after the first few time steps, high frequency oscillations of exponentially 
growing modes developed near the center of the star (i.e. at a = 0). This is most likely due to the presence of a 
coordinate singularity at the origin but also because some of the coefficients on the righthand-sides of our evolution 



system get nearly singular when moving to the center (compare the denominators of (111 with ( |T2| and ([9|). However, 
the occurrence of singular terms in general doesn't neccessarily lead to a failure of the numerical scheme. 

The solution to this problem was the utilization of additional viscous terms in the evolution equations. Since they 
do not represent any physical effect in the initial setup of the problem, they are commonly referred to as artificial 
viscosity. Here, to each of the four equations we added a Kreiss-Oliger like term of the form (see |56j for details) 

V(f)=a(D+D-+D+D~)f, (23) 

where / is a perturbation variable, a — const, is the dissipation coefficient, (a, r) are the spatial coordinates and 
D + , D~ are the standard forward- and backward-difference operators. By using dissipation coefficients oti , i = 1 ... 4 
ranging from 10 -3 — 10~ 4 it become possible to stabilize the time-evolution code against exponentially growing 
instabilities. 

As already described earlier in section [H] after a successfull simulation, the real part of the complex solution is 
taken and a discrete Fourier-transformation at several points inside the star is performed on these data to extract 
the oscillation frequencies. If N is the number of points in this time series and At is the temporal resolution, the 
corresponding frequency resolution and the Nyquist-frequency (i.e. the maximum frequency that can be resolved) are 
given by 

A/ = -^r- and f c =~ (24) 
3 NAt J 2At y ' 

Since an explicit numerical scheme has been used for time-evolution, there are certain restrictions on the absolute 
value of At. The timestep cannot be arbitrarily large and depends on the spatial resolution of our grid (CFL-criterion). 
For most of the simulations, a timestep of the order of At w 10 -6 sec and N w 10 4 was used. The total evolution 
time then is t max w 30 — 40 ms with A/ w 15 — 30 Hz and f c w 8 — 12 kHz; details are following in the next section 

EH 



IV. RESULTS 



For a spherical symmetric background and even in the slow-rotation approximation, the angular dependence of a 
mode is often described in terms of spherical harmonics Yi m . If the perturbation changes like (—1)' when applying 
the transformation r — > — r it is called an even or polar mode while odd or axial modes behave like (— 1)' +1 . In the 
axisymmetric case (i.e. for m — 0) all these modes are stable while non-axisymmetric perturbations may become 
unstable to the CFS-instability [571155] . 
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A. Axisymmetric case 

Axisymmetric perturbations, due to their simplicity have been studied in detail with perturbative methods but 
mainly via evolution of non- linear equations. 



Polar perturbations 



The results of this approach will be compared with the non-linear results published in |31j . There, a relativistic 
hydrodynamics code is used to study stellar oscillations in the Cowling approximation. The background models are 
commonly referred to as BU; they are uniformly rotating neutron stars with a polytropic equation of state with T — 2, 
K = 100 and fixed central rest-mass density p c — 1.28 x 10~ 3 in units where G = c = Mq = 1. In the nonrotating 
case, this leads to a stellar model with a gravitational mass of M = 1.4 M and a circumferential radius of R = 14.15 



km. The applied initial perturbations were decomposed into spherical harmonics and for / 
of the form 



a density perturbation 



Sp = Ap c sin 



r.(0) 



(25) 



is used as initial data. Here A is the perturbation amplitude, (r,8) denote spherical coordinates and r 8 {9) is the 
coordinate radius of the star (which is not independent of when the star is rotating) in this spherical system. In 
our simulations we mainly use a 18 x 18 spectral grid to compute the background (this already gives a very accurate 
stellar model with an accuracy of 10~ 10 ) and interpolate to our computatonal domain, which can have practically any 
arbitrary resolution; mainly we use 50 x 40, 100 x 80 or 200 x 160 grid points. 

The following figure [3] shows in the left panel a 12 msec-long section of a simulation with initial data described by 



(25 1 and a nonrotating background model on a 50 x 40 grid. The right panel depicts the logarithm of the one-sided 
power spectral density of the complete time series with a frequency resolution of A/ = 20 Hz; the data for this figure 
were extracted from the spatial position (er, r) = (0.5, 0.5). 

In the frequency plot one can see several peaks, some of them are already labelled; the vertical lines are the equivalent 
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FIG. 3: Left panel: Time-evolution of the (^-component of the perturbed 4-velocity for an axisymmetric pressure 
perturbation Right panel: The corresponding Fourier transform of the du^- time series. The fundamental radial 

mode together with a few overtones are apparent 



frequencies one can find in [31]. The strongest and one of the sharpest peaks is the one at / = 2.705 kHz and belongs 
to the fundamental quasi-radial oscillation mode (F-mode) with no nodes of the corresponding eigenfunction in 
the radial direction. Alongside with this oscillation some other modes were excited as well, most notably several 
overtones of the F-mode labelled H± , Hi and H% which have one, two or three nodes of their eigenfunctions in the 
radial direction. 

Additionally, to the estimation of the mode frequencies from a given time evolution via a Fourier transform we 
implemented a method to retrieve their corresponding eigenfunctions. The amplitude of the eigenfunction at a given 



9 



point directly correlates to the strength of the corresponding peak in the power spectral density at this particular 
point. In order to extract the eigenfunction of a specific oscillation mode one has to iterate over the computational 
domain, taking Fourier transforms at many grid points and monitor the variation in amplitude of the mode peak one 
wants to study. This eigenfunction can then be used as an improvement to the first trial eigenfunction and can be 
put back as initial data for another simulation. Usually this procedure which is called mode recycling, when applied 
repeatedly, will enhance and sharpen the mode that is recycled and will suppress additional modes that are excited. 

In figure [4] the absolute values of the eigenfunctions corresponding to the oscillation modes labelled in the previous 
figure were extracted using this technique. As already discussed earlier it is (a, r) e [0, 1] x [0, 2] and r = 1 describes 
the equatorial plane. 




FIG. 4: Top Row: The amplitude of the eigenfunctions for the F- and Hi -mode Bottom Row: The amplitude of the 

eigenfuctions for the and H3 -mode 



One can see, that 5u^ = along the equatorial plane and also that the number of radial nodes increases along the 
sequence F, Hi,H2,H 3 . 

The second strongest peak in figure [3] at / = 1.929 kHz (it is already roughly two orders of magnitude lower than 
the F-mode peak) does not belong to a quasi-radial oscillation. Instead, the extracted eigenfunction shows an angular 
dependence that is in agreement with an axisymmetric quadrupolar perturbation. By extracting the eigenfunction 
for the pressure perturbation of this mode one gets a useful initial perturbation for extracting non-radial modes and 
eigenfunctions. The result of such a simulation is depicted in the next figure [5j 

As one can easily notice, the 2 /-mode is the dominant oscillation while the first (I — 2, to = 0)-overtone, the 2 p\- 
mode, has been also considerably enhanced. In this simulation, both modes are stronger than their I = counterparts 
and several other modes which were not clearly visible in the first run (compare to figure |3| become noticable. The 
corresponding values taken from 3 1J for comparison are indicated by solid vertical lines in the power spectral density 
plot. Here we notice larger deviations than in the previous radial case but they are smaller than 5%. 

One can increase the rotation rate of the BU model and estimate how the various frequencies will change when 
the neutron star rotates faster and faster. We did this calculation for three resolutions starting from 50 x 40 and 
doubling it twice. The results for the already discussed I = -modes at the highest resolution (i.e. 200 x 160) and 
with Af — 20 Hz are summarized in Table [IT] 

the oscillation frequencies of the 2 f- and the 2 pi-mode are shown for various rotation rates. In general, 



In Table 
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these results are in good agreement with published values; their absolute differences never exceed the 5%-level. By 
increasing the resolution the most significant improvement has been observed when doubling the low 50x40 -resolution. 
For some modes the change of the frequency in these two resolutions is rather dramatic compared to the change when 
improving from the medium to the high resolution. We take this as an indication that for most purposes a resolution 
of 100 x 80 grid points suffices to get already quite accurate results. We also did convergence checks on the two 
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/(Hz) 

FIG. 5: The shape of the recovered pressure eigenfunction used for a mode recycling run and Fourier transform of 
the «£- time series with the pressure eigenfunction as initial data 



fi (kHz) 


F (kHz) 


Hi (kHz) 


H 2 (kHz) 


H 3 (kHz) 


0.0 


2.679 


4.561 


6.380 


8.178 


2.182 


2.638 


4.466 


6.255 


8.035 


3.062 


2.605 


4.435 


6.253 


8.058 


3.712 


2.570 


4.409 


6.275 


8.111 


4.229 


2.539 


4.410 


6.310 


8.156 


4.647 


2.500 


4.400 


6.330 


8.237 


4.976 


2.484 


4.392 


6.356 


8.334 


5.213 


2.456 


4.390 


6.370 


8.405 


5.344 


2.426 


4.394 


6.380 


8.411 



TABLE II: Frequencies of the axisymmetric modes F, Hi, H2 and H3 for the BU model at different rotation rates 



fi (kHz) 


2 / (kHz) 


2 Pi (kHz) 


0.0 


1.890 


4.130 


2.182 


1.890 


4.065 


3.062 


1.906 


3.970 


3.712 


1.895 


3.838 


4.229 


1.875 


3.674 


4.647 


1.844 


3.487 


4.976 


1.794 


3.275 


5.213 


1.703 


3.056 


5.344 


1.613 


2.426 



TABLE III: Frequencies of the axisymmetric quadrupolar modes 2 / and 2 pi for the BU model at different rotation 

rates 



strongest modes for I = and I = 2 at our three basic resolutions and an additional one with 75 x 60 gridpoints. The 
Iterated Crank-Nicholson scheme, which is used for time-evolution here, is first order accurate in time and second 
order accurate in space. By increasing the number of gridpoints one expects to observe a quadratic convergence in 
the perturbation functions and this was indeed the case. 

The following figure [6] shows how the frequencies of the F-, Hi-, 2 f- and 2 pi -mode change as functions of the rotation 
rate and resolution; there the line connecting circles are the values taken from [5T| . Especially the frequencies for the 
first p-mode corresponding to I — 2 agree very well with them; also the other modes show a similar behaviour. 
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FIG. 6: Top Row: The change in frequencies of the F- and Hi -mode for three different resolutions Bottom Row: 

The corresponding changes of the 2 /- and 2 pi -mode 



2. Axial perturbations 



The axial perturbations of the fluid correspond to the inertial modes, in this case small axial deviations from 
equilibrium are restored by the Coriolis force. These modes degenerate at zero frequency but have nonzero frequency 
once rotation sets in. We were able to identify some inertial modes and to compare our results with the corrsponding 
studies in |33j . In contrast to the Cowling-approximation we use here, they developed a full general relativistic 
hydrodynamics code under the assumption of a conformally flat three-metric and tested it on non-linear axisymmetric 
pulsations of rotating relativistic stars. Nevertheless, as we show, the results of the two calculations agree quite well. 

Figure [7] shows the power spectral density of the u v velocity-component of a simulation with a BU3 background 
model (in the notation of 33J ) . This neutron star has a ratio of polar coordinate radius to equatorial coordinate 
radius of r p /r e — 0.85 and rotates with an angular velocity £1 — 3.71 kHz. We chose to monitor the (^-component 
of the velocity since this is the quantity where the inertial mode signature is typically well pronounced. The three 
vertical lines in the low frequency part of the plot denote the values of the i-2, H and 12 mode frequency acco rding 
to [33]; on the right side one can identify the frequency peaks of the 2 /- and the F-mode from table pi] and 
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One can see, that there is a quite good agreement between our results and the literature values for this specific model 
and rotation rate; there is a greater deviation from the results in [33 for the lowest frequency inertial mode i-2- 
In general, the 12 mode is the strongest inertial mode in our simulation and if one compares the various r-modc 
frequencies for different angular velocities, we find the best matching results for the i\ and 12 mode and still a very 
good match for the i-2 frequencies. This is depicted in the next figure [8] 

So although we are using the Cowling approximation and not a full relativistic treatment, we can confirm the depen- 
dence of the three inertial modes i\, 12 and i-2 on the rotation rate described in the literature; a summary of our 



results is given in table IV 



B. Non-axisymmetric case 

We will now turn to the study of non-axisymmetric oscillations on rotating compact objects with emphasis on the 
m = 2-perturbations. In addition to the equation of state for the BU model series we used in the previous section 
mostly for code testing purposes, here we will apply two more equations of state. We use the polytropic parameters 
for EOS A and EOS II from [SB]; more specifically we have V = 2.46 and K = 0.00936 for EOS A and T = 2.34 and 
K = 0.0195 for EOS II. These values are given in geometric units (G = M & = 1) and with [km] as length scale. Figure 
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FIG. 7: Fourier transform of u v from a time series with BU3 as background model. The three vertical lines are from 
left to right the i-2, i\ and 12 inertial mode frequencies listed in |33] 




angular velocity / kHz angular velocity / kHz angular velocity / kHz 

FIG. 8: A comparison between the values in |33j and our results for the i\, 12 and i-2 inertial mode. Our 

simulations were performed with 200 x 160 gridpoints 



Model 


n (kHz) 


ii (kHz) 


i 2 (kHz) 


i-2 (kHz) 


BU0 


0.0 


0.0 


0.0 


0.0 


BUI 


2.182 


0.266 


0.373 


0.202 


BU2 


3.062 


0.399 


0.528 


0.285 


BU3 


3.712 


0.521 


0.651 


0.372 


BU4 


4.229 


0.602 


0.750 


0.432 


BU5 


4.647 


0.657 


0.836 


0.497 


BU6 


4.976 


0.721 


0.924 


0.539 


BU7 


5.213 


0.776 


0.971 


0.584 


BU8 


5.344 


0.790 


1.001 


0.619 


BU9 


5.361 


0.806 


1.009 


0.685 



TABLE IV: Frequencies of the three inertial modes i\, 12, i-2 for the BU model at different rotation rates 



[9] shows mass-radius-diagrams for all the EoS we are using in this paper. The black dots on every of these three 
curves are the actual models used in the simulations. The BU model is a "standard" compact object with a mass of 
M = 1.4 Mq and a circumferential radius of R = 14.15 km as one can see from the figure. In contrast to this the two 
other configurations we chose arc very close to their maximum mass limit; in particular our background model for 
EOS A has a mass of M = 1.61 M and a radius of R = 9.51 km while the EOS II model has a mass of M = 1.91 M Q 
and a radius of R = 11.68 km. They are therefore more compact and their Kepler- limit is much higher than for the 
BU model; we will see what this means for the non-axisymmetric mode frequencies in the following discussion. 
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FIG. 9: Mass-Radius relations for the EoS used to study non-axisymmetric perturbations; the black dots denote our 

actual models 



1. Polar perturbations 



The procedure for the excitation of modes is similar to the one in the axisymmetric case. Similar to the approach 
taken in [3T] we chose a I = 2- velocity perturbation of the form 



Sua = A sin 



r s (0) 



sin t) cos ( 



(26) 



with the same meaning of A and r s (9) as in ( 25 ) . Since we are working in cylindrical coordinates we have to decompose 
this perturbation into its p- and ^-component before we can insert it in our simulation. Figure [10] shows a series of 
power spectral density plots for the BU model at different rotation rates. The frontmost spectrum shows the f-mode 
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FIG. 10: Power Spectral density of the m = 2 /-mode for initial data provided by (26 1. The splitting in the 



spectrum becomes apparent for increasing angular velocities 



peak for a nonrotating star while the last one has been extracted for a star with a ratio of polar coordinate radius to 
equatorial coordinate radius of r p /r e = 0.9. In the nonrotating case the /-mode frequency is degenerated, i.e. it shows 
only one peak for the m — 2 (counter-rotating) and the m = — 2 (co-rotating) modes. This degeneracy is broken as 
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soon as rotation sets in; for £1^0 this peak splits into two peaks and breaks the symmetry between counter- and 
co-rotating modes. The very same behaviour can be observed for other modes as well. 

In the following discussion we will focus on the fundamental mode although this particular mode suffers most of 
all from the Cowling approximation (compare |31j with |33j or see |28j for an early Newtonian calculation). However 
in this paper we are interested in the evolution of the /-mode frequency with increasing angular velocity. While the 
absolute values may be incorrect by a factor of 20-30% we still should be able to make some statements about the 
characteristic behaviour of the /-mode frequency dependence on the equations of state. Figure [TI] shows our results 
for the models depicted in figure [9] In contrast to previous figures we now plot the oscillation frequencies against the 
rotation frequency which is given by v — Q/2tt. As one can see the various models have a quite different range of 




f-modc frequencies. The BU model which is the less compact also has the lowest fundamental frequency. Due to the 
model parameters the Kepler limit for this particular neutron star is reached already at vx = 853 Hz. The models 
for the other equations of state are more compact and therefore allow for higher rotation rates which can be as high 



as v = 1.758 kHz for EOS A. For each of the three EoS the m = 2 branches are those in figure 11 with the higher 
frequencies and we will see in a second how we have been able to determine this. The frequency of the fundamental 
mode scales with the compactness of the star, higher compactness means higher frequency; a property we can directly 
validate from the left panel of figure [TT] 

The right panel shows a different representation of the same picture where we plot for each model the normalized 
mode frequency tr/co against the normalized rotation frequeny v/vk where <7o is the mode frequency in the nonrotating 
limit and vk labels the Kepler limit for the rotation frequency. It is quite remarkable that although the models used 
in these simulations have very different parameters their normalized /-mode frequencies change nearly in the same 
manner when rotation is increased. It is only in the regime close to the Kepler limit where deviations for the various 
models become evident. For all rotation parameters we can write 

00 \ V K J \ V K J 



independent of the specific EoS. To determine the coefficients Ci m we made least-square fittings of the all data points 

:f 2 ] = -0.25 ± 0.02 and C$ 



we obtained from the various simulations. In particular we find = —0.25 ± 0.02 and = —0.38 ± 0.02 for the 



m = 2 solution and C { 2 % = 0.48 ± 0.03 and C { 2 % = -0.55 ± 0.04 in the m = -2 case. 

However, one should keep in mind that the results presented up to now are all derived in the corotating frame since 



this is the natural coordinate system in which our equations were formulated (see section III I . The only coordinate 
that changes when going to a stationary coordinate system is the azimuthal angle (p and the relation connecting these 
two coordinates simply is 



Vcorot — </?stat — 



(28) 
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Due to the decomposition ( 10 1 of our perturbation variables and the harmonic Fourier transformation we are per- 



forming on our numerically obtained time-series, we are effectively decomposing our time-evolution quantities like 

f r^j e iat e im f — e i(trt+m<p) ^29) 

To track a specific constant phase in time one therefore has to move on by an angle 



, corot < 

m 



after a time tQ. This means that in this case modes with a positive m are moving retrograde while waves with a 



negative m travel prograde in the comoving frame. We now insert (28 1 into (29 1 to obtain a relationship between the 



mode frequencies in the comoving and stationary frame and arrive at 

CTstat = C corot - mfl (30) 

For axisymmetric perturbations (i.e. m — 0) the two frequencies are identical; this is why we didn't start this 
discussion already in section [IV A| To track a surface of constant phase in the stationary frame a similar calculation 
like the one above yields 

stat _ (^corot ~ mft) 

<Po — z o 

m 

This means that if the frequency cr C orot is larger than mfl, then the mode is also travelling in retrograde direction 
in the stationary system. For er coro t = mVL the frequency becomes degenerate in this system while for cr coro t < mfl 
a mode travelling retrograde in the comoving frame is seen prograde in the stationary coordinate system. The next 
figure shows our results for the stationary frame. 




0.5 1 1.5 

rotation frequency / kHz 




0.2 0.4 0.6 0.8 

rotation frequency v/v K 



FIG. 12: Same as in figure [IT] but now in the stationary frame 

When we analyze the various fourier spectra in the stationary frame, we can actually see that for every of the three 
EoS the high frequency branches in figure 11 are shifted towards lower frequencies; together with equation (30 1 it 
means that this branches can be identified with the m = 2 solutions and vice versa. This is how we can identify 
the different sections of the curves. The m = 2 solutions of all models actually reach the point where a is zero; 
for the model BU this happens just at the Kepler limit, for the other models which are more compact this occurs 
even earlier. Beyond this point the /-mode is seen retrograde in the comoving frame but prograde in the stationary 
coordinate system. In this case the mode becomes CFS unstable and it can be an excellent source for gravitational 
waves. Finally, the normalized picture on the right panel now differs significantly from the corresponding picture in 
the corotating frame. This is due to the fact that the transformation from one system to the other introduces extra 
terms which are obviously model-dependent. 



1G 



2. Axial perturbations 

Non-axisymmetric axial perturbations, also known as r-modes, are known to be generically unstable to the CFS- 
instability at all rotation rates. This is an exciting class of stellar oscillations with many possible applications in 
astrophysics and gravitational wave research |13) . 

We also did a couple of simulations to specifically excite the I = 2, m = 2 inertial mode and were successful. In 
the Newtonian framework, the fundamental r-mode is of purely axial parity and thus does not mix with high order 
polar terms, see [59] • It also has been shown in [BO], that the contribution of higher order terms introduced by the 
fully general relativistic treatment can be neglected in the case of slow rotation. As for the polar non-axisymmetric 
oscillations we start with a perturbation of the ^-component of the 4- velocity and expect a dependence ~ 1/r sm(6) 
of Sug, see also |61j . Typically, several mode recycling runs are needed to get a sharp and clean signal in the power 
spectral density. Also, due to the low frequencies of inertial modes especially at small rotation rates one needs much 
longer evolution times than for pressure driven modes. We chose to cancel the time-evolution after roughly 50 ms to 
extract the eigenfrequencies. The numerical code is still stable there and can in principle evolve the initial data for a 
longer time intervall, leading to a more accurate frequency determination. However, we found that an evolution time 
of about 50 ms and a spatial resolution of 200 x 160 gridpoints is already quite good for a first estimation; there is 
only a marginally change when using longer integration intervals. We compare our code with results for a BU6 star 
rotating at 93% of its maximum speed as described in [62 where a nonlinear general-relativistic code has been used 
to study the saturation amplitude of r-modes. The following figure [13] gives a summary of the results. 
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FIG. 13: Top Left: Time evolution of Su p Top Right: Corresponding power spectral density plot Bottom Left: Shape 
of the Su p eigcnfunction Bottom Right: Shape of the Sur cigenfunction 



By means of repeated mode recycling, the fundamental r-modc can be significantly enhanced as depicted in the 
PSD plot of figure [13] We also checked the angular dependence of the extracted eigenfunction. Since 5ug should be 
~ sin(#), the p- and ^-components of the perturbed 4- velocity are modified by an additional factor of cos(0) and sm(6) 
respectively. Keep in mind that the grid variable r which takes values from 1 to 2 can be thought of similar to the polar 
angle 9 in spherical coordinates; see the discussion of figure [T] in this paper. Then it is clear that the eigenfunction for 



5u p and Su^ depicted in the bottom row of figure 13 show indeed the expected behaviour. The fundamental r-mode 
has a frequency of f c — 518 Hz in the comoving frame which translates to a frequency of /, = 1.066 kHz in the inertial 
frame. This is in excellent agreement with |62j where they found it at /, = 1.03 kHz and also with |63j where they 
saw the mode at fa = 1.05 kHz; nonlinear effects in the simulations of [B2] may explain the larger discrepancy there. 
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V. CONCLUSIONS 



In this work we present a study of the oscillation properties of fast rotating relativistic stars for bot both axisym- 
metric and non-axisymmetric perturbations. The study was based on the Cowling approximation using a 2D version 
of the perturbation equations. The results for axisymmetric perturbations are in excellent agreement with earlier ones 
while the non-axisymmetric results are the first of their kind in the literature. We demonstrated the neutral points 
for the onset of the CFS instability and suggested possible normalizations which can be used in order to extract the 
parameters of the rotating star. 

The method presented here (evolution of the 2D linearized equations) can be extended by including the perturbations 
of the spacetime. This will offer the possibility in testing the earlier results [B3] for the onset of the secular instability 
in fast rotating stars. Moreover, the dependence of oscillations frequencies on rotation will be based on the exact 
results and will not rely on the Cowling approximation. Finally, the effect of differential rotation on the spectra can 
be studied both for testing earlier, mainly Newtonian, results [43, 65 but more importantly in finding the correct 
dependence of the frequencies on rotation as well as the neutral points for the onset of the secular instabilities. This 
is actually an important step since newly born neutron stars are expected to rotate differentially at least during the 
stages that they will be in an oscillatory phase or the even when they are secularly unstable. 
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